function VBas = affine_vbasis(xv, yv, names)
%AFFINE_VBASIS Velocity field bases of affine flows
%
%   The entire Affine flow family has six base flows, whose Lie algebraic
%   representation are respectively given below:
%
%   * x-translation (xt): [0 0 1; 0 0 0];
%   * y-translation (yt): [0 0 0; 0 0 1];
%   * x-scale (xs):       [1 0 0; 0 0 0];
%   * y-scale (ys):       [0 0 0; 0 1 0];
%   * rotation (ro):      [0 -0.5 0; 0.5 0 0];
%   * shearing (sh):      [0 0.5 0; 0.5 0 0]; 
%
%
%   VBas = AFFINE_VBASIS(h, w);
%
%       Generates the entire set of affine velocity field basis, of 
%       specified size.
%
%   VBas = AFFINE_VBASIS(xv, yv);
%
%       Generates the entire set of affine velocity field basis, with
%       provided coordinate vectors.
%
%       AFFINE_VBASIS(h, w) is equivalent to AFFINE_VBASIS(1:h, 1:w);      
%
%   VBas = AFFINE_VBASIS(h, w, names);
%   VBas = AFFINE_VBASIS(xv, yv, names);
%
%       Generates a subset of bases (specified by short-names).
%
%       For example, AFFINE_VBASIS('xt', 'yt', 'ro') generates three
%       bases, respectively for translation and rotation.
%
%       AFFINE_VBASIS() is equivalent to 
%       AFFINE_VBASIS('xt', 'yt', 'xs', 'ys', 'ro', 'sh');
%

% Created by Dahua Lin, on April 5, 2012
%

%% verify input arguments

if isscalar(xv)
    xv = 1 : xv;
else
    if ~(isfloat(xv) && isvector(xv))
        error('affine_vbasis:invalidarg', 'The 1st argument is invalid.');
    end
end

if isscalar(yv)
    yv = 1 : yv;
else
    if ~(isfloat(yv) && isvector(yv))
        error('affine_vbasis:invalidarg', 'The 2nd argument is invalid.');
    end
end

if nargin < 3
    names = {'xt', 'yt', 'xs', 'ys', 'ro', 'sh'};
else
    if ~iscellstr(names)
        error('affine_vbasis:invalidarg', ...
            'The names should be a cell array of strings.');
    end
end

%% main

K = numel(names);
[X, Y] = meshgrid(xv, yv);
[h, w] = size(X);

coords = [X(:).'; Y(:).'];

Bx = zeros(h*w, K);
By = zeros(h*w, K);

for k = 1 : K
    
    [Y, u] = get_liealg(names{k});
    
    vv = bsxfun(@plus, Y * coords, u);    
    vx = vv(1,:).';
    vy = vv(2,:).';
    
    Bx(:,k) = vx;
    By(:,k) = vy;
end

VBas.tag = 'vbasis';
VBas.size = [h, w];
VBas.K = K;
VBas.Bx = Bx;
VBas.By = By;


%% Lie algebra matrices

function [Y, u] = get_liealg(name)

switch name
    case 'xt'
        Y = zeros(2,2); u = [1; 0];
    case 'yt'
        Y = zeros(2,2); u = [0; 1];
    case 'xs'
        Y = [1 0; 0 0]; u = [0; 0];
    case 'ys'
        Y = [0 0; 0 1]; u = [0; 0];
    case 'ro'
        Y = [0 -0.5; 0.5 0]; u = [0; 0];
    case 'sh'
        Y = [0 0.5; 0.5 0]; u = [0; 0];
end

